Cosmic propagators at two-loop order 
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We explore the properties of two-point cosmic propagators when Perturbation Theory (PT) loop 
corrections are consistently taken into account. We show in particular how the interpolation scheme 
proposed in [T] can be explicitly used up to two-loop order introducing the notion of regular parts 
for the contributing terms. Extending the one-loop results, we then derive and give semi analytical 
forms of the two-loop contributions for both the cosmic density and velocity propagators. These 
results are tested against numerical simulations and shown to significantly improve upon one-loop 
results for redshifts above 0.5. 

We found however that at lower redshift two-loop order corrections are too large partly due to a 
strong sensitivity of those terms to the small scale modes. We show that this dependence is expected 
to be even larger for higher order loop corrections both from theoretical investigations and numerical 
tests, the latter obtained with Monte Carlo evaluations of the three-loop contributions. This makes 
small-scale regularization schemes necessary for the use of such higher order corrections. 

PACS numbers: 98.80.-k, 98.65.-r 



I. INTRODUCTION 

Precision measurements of the statistical properties of the large-scale structure of the universe is the key to current 
large projects of observational cosmology. They aim in particular at measuring the galaxy power spectrum at the scales 
where Baryonic Acoustic Oscillations can be detected ([2HS]) or at measuring the projected density power spectrum 
from weak lensing observations (as stressed in [B] and which is the part of the core program of the Euclid mission [7]). 
In this context, semi-analytic tools are now developed that aim at computing, from first principles, the statistical 
properties of the cosmic fields at large-scale. Standard Perturbation Theory approaches (see [5]) are however not 
very efficient in producing well behaved next-to-leading order terms to the spectra or bispectra. It leads indeed to a 
resummation series that has poor converging properties (see [9] for more insights into that). Alternative approaches 
have been recently developed that propose, implicitly or explicitly, a reorganization of the perturbation series. This 
is the case in particular for the Renormalization Perturbation Theory (RPT) approaches developed initially in [§]. 
Other approaches are based on the writing of closed system of integro-differential equations as in the closure theory 
approach, [ID] , or in the Time Flow equation hierarchy presented in [TT]. We put ourselves in the context of the RPT 
approaches and its extensions which are based on the introduction of the two-point propagator and on a re- writing of 
the series expansion with the help of that object. A related but alternative approach has been proposed in [12 where 
was introduced the concepts of multi-point propagators. Interaction vertices can then be formally renormalized as 
well. It has been shown there that a new scheme for the computation of spectra and poly-spectra can be obtained 
form the multi-point propagators, from the the so-called V expansion. 

The focus of this paper is precisely on the behavior of the two-point propagator and on the precision one can 
theoretically attain in its description. In particular one wishes to evaluate the importance of the multi-loop corrections 
to its expression, in practice up to second order, and to see whether these corrections can be safely computed. 
Comparisons with N-body codes will be used taking advantage of the interpolation scheme proposed in pQ. In this 
paper it is indeed shown how multi-loop corrections can be consistently used to get predictions on multi-propagator 
expressions with theoretically increasing precision. These investigations aim in particular at pinning down the validity 
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range of this type of calculations. As we will see in the following the two-point propagator is the quantity, at a given 
order in PT, that is the most sensitive to high-fc modes in the computation of power spectra. Quantifying the 
importance of the small scale physics, for large-scale predictions, is a daunting quest. Some phenomcnological results 
have been proposed in [T31 Q3] to investigate for instance the importance of shell crossings. This is also the idea behind 
the works in [T3 [TS] where the impact of the small scale physics is tentatively introduced as extra source terms. In 
this paper we eventually investigate this problem from a PT approach. 

The paper is organized as follows. In Sec. [IT] we recall the motion equations for gravitational instabilities of a cosmic 
fluid. In Sec. |III| we compute the two-loop corrective terms and propose analytic fits for them. In Sect. |IV| we first 
estimate the contributions of the three-loop terms from Monte-Carlo computations. Our results are tested with the 
help of consistency relation we provided in Sect. |III| We then explore in more detail the expected properties of the 
loop corrections at arbitrary orders introducing the notion of kernels and henceforth quantifying the sensitivity of the 
results to large and small scales modes. Our concluding remarks are presented in Sec. [V] 



II. THE COSMIC TWO-POINT PROPAGATORS 



A. The equations of motion 



We are interested here in the development of cosmological instabilities in a cosmological dust fluid. In general 
the dynamical evolution of such a fluid can be described with the Vlasov equation. As usual we then restrict our 
investigations to the regime where multi-flow regions play a negligible role but we will discuss this assumption later 
on. In the one flow limit, the motion equation then takes the form of a set of three coupled equations relating the 
local density contrast, the peculiar velocity field and the gravitational potential (see [8]). 

At linear order these equations can easily be solved for an arbitrary background cosmology. One generically finds 
a growing solution and a decaying solution. Let us denotes D + (r) the growing mode solution for the density contrast 
/+( T ) its logarithmic derivative with respect to the expansion so that, 

S(k,f ] ) = D + (n)8 (k), 6(k, V )/H = -f + ( V )D+(rj)5o(k) (1) 

is the solution for the growing mode and similarly, 

S(k, v ) = D.(r,)S Q (k), 9(k, V )/n = -/_fo)D_(»j)<Jo(k) (2) 

for the decaying. 

Following [17], the motion equations describing a pressure-less fluid in the one-flow limit can be written in a compact 
form with the use of the doublet \& a (k, r), defined as, 

¥ a (k,r) = (5(k,r), -_L_0(k,r)), (3) 

where H = dhia/dr is the conformal expansion rate with a(r) the cosmological scale factor and r the conformal time 
and where the index a = 1, 2 selects the density or velocity components and which makes explicit use of the growing 
solution. 

It is then convenient to re-express the time dependence in terms of the growing solution and in the following we 
will use the time variable rj defined as 

r) = log D+(ri) (4) 

assuming the growing factor is set to unity at initial time. Then the fully nonlinear equations of motion in Fourier 
space read [3] (we henceforth use the convention that repeated Fourier arguments are integrated over and the Einstein 
convention on repeated indices), 

d 

— * a (k,r/) +O ab (r/)* fc (k,77) = 7 abc (k, ki, k 2 ) * 6 (ki, 17) * c (k 2) rf), (5) 



where 



-1 

3 n m 3 n m 



(6) 
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and the symmetrized vertex matrix j abc describes the non linear interactions between different Fourier modes. Its 
components are given by 

. , s c , . x |k X + k 2 | 2 (k! • k 2 ) 

722 2 (k,ki,k 2 ) = d D (k-ki - k 2 ) - 



), and 7 = otherwise, where 5v denotes the Dirac delta distribution. The matrix j abc 
is independent on time (and on the background evolution) and encodes all the non-linear couplings of the system. 
The formal integral solution to Eq. ^ is given by (see [9l \17\ [18] for a detailed derivation) 



* Q (k,7j) = g ab (r,) &(k) + / d»/ff o6 (»7 > t/)7W,(k > k 1 ,k 2 )* c (ki 1 »/)* ll (k ai »/) 1 (8) 

Jo 

where </> a (k) = ^ a (k,ri = 0) denotes the initial conditions, set when the growth factor D + = 1 and where g a b{ii) is 
the linear propagator that is Green's function of the linearized version of Eq. ([5| and describes the standard linear 
evolution of the density and velocity fields from their initial state. 

In the following calculations we will be using the value of the Vl ab matrix to be that of the Einstein de Sitter 

3/2 

background that is effectively assuming that D_ scales like D + . This is known to be a very good approximation 
even in the context of a A— CDM universe (see for instance [15] for an explicit investigation of the consequences of 
such an approximation). 

The ensemble average of any quantity can then be built out of those of the initial fields. They are entirely defined 
from the initial power spectrum of density fluctuations P a b{k) define as, 

<0 a (k)0 b (k')) = <5 D (k + k')P a6 (fc). (9) 

In what follows most of the calculations and applications will be made assuming initial conditions in the growing 
mode, for which </> a (k) = (5 (k)w a with u a — (1, 1), and therefore with P ab (k) = Po(k)u a u b . 



B. Definition and general properties of the two-point and multi-point propagators 

The g a b(s) appears clearly to be one of the key ingredient of these constructions. It can be described as the the 
linear propagator giving the variation of the mode amplitude as time is running. The idea at the heart of the RPT 
approach is to generalize this operator beyond linear theory. More specifically the quantity ^^fej expresses the 
way Vl/a(k, s) depends on </>b(k') as a function of time s. This function however depends on the stochastic properties 
of the fields. One can nonetheless define its ensemble average, G ab (\t), as 

2^) = aD( k- k ')^( M ). (10) 

This quantity depends on the initial fluctuations through the mode couplings. The ensemble average is made precisely 
over these modes. The Dirac-(5 function it leads to is due, as usual, to the homogeneity of the underlying statistical 
process. More generally one can define the {p + l)-point propagators T^ pS> as (see [H]), 

? ( »,":.Xw ) = Mk - r °u* ("> 

The Dirac ^-function that appears in this expression is a direct consequence of the assumed statistical isotropy of 
space. Note that the propagators can be defined for any time interval. We will however use it here assuming the early 
time to be the initial time for which one can assume that the cosmic modes correspond to the Matter Dominated era 
growing mode (see |20j for explicit investigations of the validity of such an approximation) . In computing the power 
spectrum or bispectrum, it is convenient to introduce the reduced propagators as 

r^Od, . . . , k p ; V ) = rit.^ki, . • • , k p ; V )u bl ...u bp . (12) 

Our main interest here is the late time evolution of the propagators keeping only the most growing modes contri- 
butions. In particular, in the following we compute the expression of G ab order by order in perturbation theory, up 
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to two loop order. Such results can be obtained from a formal expansion of ^/ a (k, 77) with respect to the initial field, 

00 

* a (M)=x;*i n) (M) as) 



71 = 1 

with 



¥<T>(M) = /" d 3 k 1 ...dXfc(k-ki...n)^^... 6n (ki,...,k«;»7)06 1 (ki)...K(kn) (14) 

where J-^ are fully symmetric functions of the wave-vectors. Note that these functions have in general a non-trivial 
time dependence because they also include sub-leading terms in 77. Their fastest growing term is of course given by 
the well known {F n , G„} kernels in PT (assuming growing mode initial conditions), 

■ 7 o6i6 a ...6„( k i' • • ■ ,k„;?7)u 6l ...u bn = exp(ray) i^ n) (ki, ..,k„) 

with Fi — F n and F 2 = G n (density or velocity divergence fields respectively) . Those quantities can obviously be 
computed in a standard Perturbation Theory approach. Formally one can then writes, 

p(p) - r (p) +r (p) +r (p) +r (p) + 

1 a — 1 a, tree ~ 1 o,l-loop T 1 a,2-loop T 1 a,3-loop 

When keeping, for each term, the only fastest growing mode, the expressions for the perturbative corrections, „_ loop 
can be rather simplified and written in terms of the standard PT kernels. Those expressions can be written more 
specifically for the two-point propagators. Let us define then the 2 propagators Gi+ and G 2 + as, 

Gi+(M) = G la (k,r))u a = G n (k,r)) + G 12 (k,ri) (16) 

G 2 +(fc,??) = G 2a (k,v)u a = Gai(M) + G 22 (fc,?7), (17) 

for which we have, 

G*+ loop (fc, 77) = 3e 3 "|d 3 qFi 3 )(q,-q,k)P (<z) (18) 

G?; loop (*,»7) = 15e 5 " y'd 3 q 1 d 3 q 2 ^ 5 )(q 1 ,-q 1 ,q 2 ,-q 2 ,k)P ( (?1 )Pofe) (19) 

with qi = |qj| and where the function F^ 1 are the symmetrized PT kernels of the ra-th order. One can generalize 
the above expressions to the one at a general n-loop order [35], 

G^ loop = sip e^+^ J d 3 q 1 ...d 3 q, i Fi 2 "+ 1 )(q 1 ,-q 1 ,...,q„,-q n ,k)F (< Zl )...F ( (Z „) (20) 



with a symmetry factor given by 



3 W = ^n^( 2 2 M=(2n + l)!!. (21) 



Note that the explicit expressions for the standard PT kernels can be obtained from the recursion relations (see also 
0), 

Fi^kx) = {1,1}, (22) 

(n) 7i—l 

Fi"'(ki,...,k„) = ^L. ^ 7abc (q 1 ,q 2 )^(™)(k 1 ,...,k m )^"- m) (k n _ m+1 ,...,k„)+sym. (23) 



m—l 



(n) 

with qi = kj. + ■ • • + k m and q 2 = k m+1 + • • • + k„. Here the matrix a ab is given by 



<>) = 1 f 2 ™ + 1 2 ) (24) 

ab - (2n + 3)(n-l) \ 3 2n ' ' U4) 
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where all the other terms obtained by permutations of the wave modes are included (there are n! of those.). 

The resulting expressions depend on the actual linear power spectra. In this sense the result is a priori model 
dependent. Note that in practice we also assume that the structure of the coupling kernels are that of an Einstein de 



Sitter universe. In general the coefficients appearing in the expression of Eq. (24) are time dependent. In practice 
they appear to change very weakly. As mentioned before, that does not restrict much the scope of these calculations. 

The propagators are of further interest because of their high-/c behavior that can be derived explicitly taking into 
account large-scale modes in a full nonlinear manner. Those properties were first derived in [SJ 121) from diagram 
resummations. They were further engrained in a more general framework in |22j with the introduction of the so-called 
eikonal approximation we will frequently refer to in the following. From these calculations we expect to have in the 
high-/c limit, 

G ab (k, V ) -to^exp (-^(e" - l) 2 ) , (25) 

where cxd is the r.m.s. of the displacement field |37j . 

9 f d 3 k , N 

°3 = J w Po(k) - (26) 

This expression however does not reproduce the one-loop expression of the propagators. Let us denote SG^T oop the 
expression of the one-loop correction to G a b and 6G 2 ^ loop its expression at two-loop order. As shown in [9] , it is possible 



to build an interpolating expression for G a b that reproduces both the asymptotic property ( 25 ) and the one- loop order 
result and which is valid at any time and for any wavenumber k. The construction of the interpolation function is 
however not systematic and relies on specific properties, k and time dependences, of the one-loop expressions. An 
alternative form for constructing a systematic interpolation functions of multi-point propagators has been proposed 
in pp. It makes then possible to incorporate corrective terms of arbitrary order into the expression of the propagators. 
We will make use of these results in the following. 



C. The one-loop order results 



The expressions of the propagator have been computed in |21j at one-loop order. They appear to be a sum of four 
terms each of which can be factorized in time and k dependence. Here and in the following we are interested in the 
most growing part of the propagators only and for initial growing modes. In this regime, we then have, 



Gi; loop (M) = 47re 3 " q 2 dq P (q) f(k,q) 



G 



1 — loop i 

2+ 



(k, V ) = 47re 3 " / q 2 dqP (q)g(k,q) 



with 



f(k,q) = 



g(k,q) = 



3(2fc 2 + 7 g 2 ) (fe 2 -g 2 ) 3 log 



(k-qf 
(k+qf 



4 (6k 7 q - 79k 5 q 3 + 50k 3 q 5 - 21kq 7 ) 



3(2k 2 +q 2 ) (k 2 ~q 2 f log 



(k-q) 2 
(k+qp 



2016fc 3 <? 5 

4 (6k 7 q - 41k 5 q 3 + 2k 3 q 5 - 3kq 7 



672fcV 



(27) 
(28) 

(29) 
(30) 



The functions f(k, q) and g(k, q) are dimensionlcss functions (e.g. they depend only of k/q) that are entirely determined 
by the vertices and by the linear propagator. More precisely we have, 



f(k, Q) = ~J d 2 <? F 3 (k, q, -q), g(k, q) = i- J d 2 q G 3 (k, q, -q), 



(31) 



where q is the relative angle between q and k. 

A few remarkable properties are worth mentioning. In particular the low q behaviors of f(k,q) and g(k,q) are 
directly related to the exponential cut-off of large-fc. We then have 



f(k,q) ~ g(k,q) 



lfc 2 
6^ 



when q <C k. 



(32) 



(> 



and the same property is true for g(k,q). The low k properties of f(k,q) and g(k,q) are less universal. There is 
however an important result which is due to the low k behavior of the vertices. It implies that f{k,q) ~ k 2 /q 2 and 
g(k,q) ~ k 2 /q 2 when fc<g. These latter properties ensures a rapid convergence of the integrals in Eqs. ( 27|28 1 and 
are thus important form a practical point of view. More precisely we have 



In the next Section we will see how these results are changed for the two-loop results. 



D. Multi-loop propagators: constructing regular parts 

Clearly propagators at each order are expected to exhibit specific behavior at large wave number. The leading k 
behavior is thus expected to behave like, 

i (_l.\2p 

GST^M) - ^SNdV - i) 2p - (34) 



This is precisely what has been obtained in [9] and allows the construction of the form ( |25| ). Furthermore, as shown in 
P] it is possible to write interpolating forms that incorporate both the standard PT results and the expected large-/c 
behaviors of the propagators. At two-loop order this form for the propagator reads 

G afc (fc, r,) « (g ab (r,) + G^ loop (fc, 77) + ■ ■ ■ + G^ loop (k, „) + c.t.) exp (-^(e" - l) 2 ) (35) 

where c.t. is a counter term. It is such that the p-loop order of this form reproduces the exact expression of G a b at 
the same order. For instance a simple calculation leads to 

ct. = 9ab (v)^f^ if + \gM (^fj (e" - If + " l) 2 G^ Ioop (fc, r?). (36) 

for the expression of the counter term up to two-loop order (the first term being its expression at one-loop order). 

Here we would like to take advantage of the general properties of the propagators, in particular on its dependence 
with <7d, to introduce consistency relations that can be used in practice to test Monte-Carlo results. The idea is to 
see G p ^ loop (k, 77) as a functional of the linear power spectrum P{k) and to express the dependence of G p ^ loop (fc, 77) 
with respect to <j&. Let us do this construction order by order. We can define the regular part of G^ loop (k, rj) by 

G^°° P (M) = -y^(e" - lf 9ab (v) + Res G^ l0 ° P (M). (37) 

In this expression Reg G* b lo ° p (fc, 77) is the part of G^ loop (fc, 77) which remains finite even when one lets <7d go to infi nity . 
The next step is to define the regular part of the two-loop results. The <j\ term has already been given in Eq.(34|. 

The <j\ term corresponds to the one-loop correction of Rcg G* fc l °° p (k, 77) taken in the high-fc limit, 

G^ lo ° P (M) = - l) 4 ft*fo) ~ yofc" " l) 2 Rog G:- lo ° P (fc,77) + R -G^ lo ° P (fc,77), (38) 

where again Rog G 2 b loop (fc, 77) is finite in the formal limit cr^ — > 00. As we will explicitly use it in the following we go 
one step further for the three-loop term, 

G^ lo ° P (k,v) = - lf 9a M + ^a d (e" - l) 4 R ^ lo ° P (k, 77) 

,2 (39) 
-y^(e" If ^ l0 ° P (k,v) + Rcs G^°° P (M). 

In this expression the first term is obtained from the large k behavior of the three-loop expression, the second from 
the large-fc behavior of the four-point propagator and the third from the six-point propagator computed taken in 
configurations where all modes are in the high k limit. 
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The expression (39 1 can be used as a consistency checks for Monte Carlo integrations of three-loop contributions. 
Indeed one expects the last term Rog G^ b loop (fc, 77) to then remain finite. And from a a& expansion point of view the 



expression (35) simply reads, 
G 



ab (k,v) « ( 9 aM + ^G^M) + ^G 2 ab lo ° P (k, v )) exp (-^(e" - l) 2 ) . 
This form underlines the importance of the regular part of each contributing term we have defined. 



(40) 



Note that in the previous calculations the expression of <7d is computed a priori from linear theory, i.e. Eq. (26), but 



its computation can incorporate in principle higher order terms. This is a direct consequence of the eikonal results 
where it is shown that the variance of the displacement field can be computed beyond the linear regime [38 • For 
instance the 1-loop correction to the displacement field can be explicitly computed and it reads, 



(4*) 



dq 1 



.1.1-1...,. - ll2s , l(j / -^-PM^-PM 



{qi-ql 4 log 



(gl + 92) : 
(Qi ~ Q2) ' 



MlQ2 (qf + q 2 2 ) (qf - + qf) 



(41) 



In such a case the corrective term can be taken into account in the expression of the regular part. Eq. (37) is then 
left unchanged but Eq. ( 38 ) is modified into 



Gl? oop (k,v) = -* r 4x-ioope 2 'V - ifgaM + ^-< lin .(e" - ifgM - ^< lin .( e " - if R ^G^ lo ° P (fc, v ) 



Re gr . 2 - lo °P 



(M)> 



(42) 



in order to incorporate the one- loop correction that appears in the counter term of Rcs G* fc loop (A:, rj). Then Eq. (40) 
is changed into, 



(g ab (v) + Rcs G 1 ab loop 



(k, V) + Rcg G:- lo ° P (fc, „)) exp (- fc2(gili "- + ^ 1 - loope2,?) ( e " - 1) 



(43) 



when the propagator is written up to two-loop order. Note that the extension to the three-loop terms is a bit 
more cumbersome. This is due to the fact that it is then necessary to take into account fourth order cumulant, the 
trispectrum, of the displacement field that appear at this order of PT calculations which should also be incorporated 
in the exponential fact or [3^] as explicitly shown in |23j . 



The importance of the change from ( 40 ) to ( 43 ) will be discussed in the next section from the explicit computations 



of the two-order expression. As we will see the changes it introduces are very benign. In practice we use an exponential 
factor which is computed either from linear theory or from explicit measurements in simulations. It is also to be noted 
that the concept of regular parts can be extended to other types of PT diagrams, not necessarily those contributing 
to the two-point propagators. We leave this analysis for a further study. 



III. TWO-LOOP RESULTS 



A. The contributing diagrams 



-O 



FIG. 1: Diagrammatic representation of the linear propagator. <j>b represents the initial conditions and g a b is the time dependent 
propagator. This diagram value is given by Eq. Q when the second term of the right hand side has been dropped. 



The computation of G 2 ab loop is still an arduous task that relies on the computation of many contributing terms. 

(5) 

These terms can be formally obtained from the expression of that is the expression of the cosmic fields at fifth 







FIG. 2: Diagrammatic representation of the fields at second order. This diagram value is given by Eq. |8f when one replaces 
\P C and 'I'd in the second term of the right hand by their linear expressions. In the diagram, each time one encounters a vertex, 
a time integration and a Dirac delta function in the wave modes is implicitly assumed. 



Y' 5 >(k) = 



FIG. 3: Diagrammatic representation of the fields at fifth order. Three different diagrams are found to contribute. 

order in the initial conditions. Formally the resulting expression G^ loop can be written, 

G?; loop (M) = (4tt) 2 e 5 " J ql& qi P { qi ) J q 2 2 dq 2 P (q 2 ) q\,q 2 ) (44) 
G 2 2 - lo ° P (k, V ) = (47T) 2 e 5 " J q 2 d qi P ( qi ) J q 2 2 dq 2 P (q 2 ) G(k, qi ,q 2 ). (45) 

where the expressions of the functions J- and Q can be derived from F§ and G5 respectively, 



(47T) 



X — J d 2 qi d 2 q2 F 5 (k,q 1 ,-q I ,q2,-q 2 ) > (46) 



g{k, qi ,q 2 ) = J d 2 qi d 2 q 2 G 5 (k,qi,-qi,q 2 ,-q 2 ), (47) 

Interestingly these expressions can be represented at a diagrammatic level. 

In Fig. [T] defines the linear propagator. In Fig. [2] we show how the coupling term is represented. In this diagram 
we show the formal expression of the second order fields. And Finally, in Fig. [3] we represent the 3 types of diagrams 
that contribute to the fifth order expression of the fields. It is this one that one should use to construct the two-loop 
contributing terms. They are obtained when 2 of the incoming lines are "glued" together: each time it introduces a 
linear power spectrum factor represented by a crossed circle. We are then led to the contributing diagrams depicted 
on Fig. [1] At this stage diagrams represent a mere transcription of the expansion of the motion equations. A detailed 
description of the procedure to draw the diagrams and compute their values can be found in [9] . 

There appear to be 9 different diagrams contributing to J- and Q. Note that each of these diagrams comes with a 
symmetry factor. The computation of these symmetry factors can be done by counting the number of possibilities one 
has to get the same final expression out of a given contributing diagram to 'J/' 5 ). Finally the symmetry factors of the 
diagrams of Fig. [3] are respectively, 2, 1, 2, 4, 1, 1, 1, 1, 1 for diagrams 1 to 9. Note that the diagrams contributing to 



the one- loop correction of the displacement field, Eq. (41 1, are the diagrams 1, 2 and 5 (with their previous symmetry 
factors) when computed in the eikonal limit. 



A priori the 2 angular integrations contained in Eqs. ( 46|[47 l can be done irrespectivevely of the details of the shape 



of the power spectrum and can therefore be done in a model independent way. The whole computation is however 
cumbersome and leads to a very large number of terms, the complexity of which depends on the diagrams. Some 
diagrams correspond to "simple" diagram (the diagrams 1, 5, 7 and 9) the expression of which can be expressed in 
terms of the function / or g only. This is not the case however for the others. Such angular integrations involve 
various sets of functions. One is built out of 



and another of 



d»j 1 = 7Tl0 g (-W + /W-) 

y |k x + q| 2 |k 2 - q| 2 k 3 q ^(k 2 2 - q 2 ) k\ + q 2 {-k 2 + k 2 + q 2 ) 



9 




with 

W± = ±2(kj - q 2 ){k 2 - q 2 ) ± 4fc|g 2 + 4k 3 qyjkjk% + g 4 + q 2 (k 2 - k\ - fc2 2 ). (50) 
These two functions are those that appear in the expression of the one-loop correction to I^ 2 '. The calculation of T 
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and Q involve however further integrals of the form, 



(for diagram 5), 



(for diagram 6), 



(for diagrams 7 and 8) or 



d 2 gid 2 g 2 i rm rm rrr (51) 

| qi +k| 2 |q2 + k| 2 | qi + q 2 + k| 2 



d & d & I T 141 T TT^2 ( 52 ) 

qi + q2| 4 |qi + q2 + k| 2 



& 2 qi& 2 q 2 | ■ pr. . ; TT7? ' (^3) 

|qi + qal |qi + k l lqi + q2 + k| 2 



d2| ?id 2 <? 2 i TTl2i Ti2i , 12 ( 54 ) 

|qi +k| 2 |q 2 -k| 2 |qi +q 2 | 2 



(for diagram 9). The angular integrations in those expressions can be done only partially using the relations (48 49 1. 
We are left with expression that involve one single angle (that could indifferently be the one between q! and q 2 or k 
or the angle between q 2 and k depending on the choice of order of angle integrations) . The resulting expressions are 
yet far too complicated to be reproduced here. In the following we rather give some general asymptotic properties 
and propose some fitting functions. 



B. The asymptotic behaviors 

The high k asymptotic properties of T and Q can be inferred from the fact that G^ loop should be consistent with 
the second order term of expression (251 in a perturbative expansion. That implies in particular that 



This asymptotic property can actually be extended in a more subtle way. Let us assume that both k and q 2 are finite 
and that q\ is much smaller than both q 2 and k. In this case T(qi, q 2 , k) is nothing but the one-loop correction of I^ 3 ) 
taken in a particular configuration. In this limit though the one-loop correction is independent on the configuration 
and is given by — k 2 /(6q 2 ) times the tree order result. As a consequences one gets a relation between T and /, 

k 2 

T(k,qi,q 2 ) -> -— ^f(k,q 2 ) when q x < k and q x < q 2 . (56) 

A similar expression naturally holds when the roles of q\ and q 2 are inverted. These two results can be recapped in 
the following form, 

J r (fc,qi,g 2 ) -> 7}f(k, qi)f(k,q 2 ) when q x < k or q 2 < k. (57) 

Note that The validity of this approximation is controlled by the magnitude of f(k,q±) + f(k,q 2 ). Moreover when 
either q x or q 2 is large one expects to obtain either a 1/q 2 or l/q 2 behavior respectively. This comes from the IR 
properties of the kernels. All these properties can be encapsulated in the following forms, 

1 k 2 
F(k,qi,q 2 ) = xf(k,q x )f(k,q 2 )- 2 2 a f {q 1 /k 1 q 2 /k), (58) 
* 9i + 1 2 



1 k 2 

G(k,qi,q 2 ) = ^g{k,qi)g{k,q 2 ) 2 2 a g (q 1 /k,q 2 /k), (59) 

1 Qi + Q 2 



where ctf and a g are both finite. The first term dominates when either q x or q 2 is small, the second when both are 
large. Note that the second term is always negative. This is illustrated on Fig. [5j 
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The functions a/ and a g denned in Eqs. (58 59) are shown on Fig. km They are clearly finite in all possible q\ and 
q 2 ranges. It is actually to compute their asymptotic expressions whenlc is small. They read, 

a/ (9i,32) Af(<7i,92) and a g (q 1 ,q 2 ) -> j3 g (qi,q 2 ) (60) 

Pf (.01,92) = 27518400g| g 7 {^2 (5760 9l 10 + 19365g 2 g 8 - 114653^ - 114653^ + 19365g 2 8 ^ + 5760g : 
+15 (384^ + 2699^? + 384g 4 ) {q\ - q 2 2 ) 4 log [fc^T 



with 



} 



(61) 



and 

Pg (?i,?2) 



Qni£»tml r l 4 ^ (I7280gi° + 120495«§g? - 572759g 4 ^ _ 5 72759 9 fg 4 + 120495g 2 8 9 2 + 17280?. 
302702400<7{ q 2 L v 

r (gi - 92) 2 
(9i +92) 2 



ION 



ouz ( vz,<ivuq 1 q 2 y. 

+ 15 (I152g 4 + 12257 9 ^g 2 + 1152<? 4 ) 



9? -9l) 4 log 



These forms are reached whenever both gi/fc and 52/^ are large 



(62) 
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FIG. 6: The functions af(k, qi, q2) and a g (k, qi, 52). Each plot shows the superposition of the result of the numerical integration 
and the fitted functions. The results of the numerical integrations show numerical instabilities that the fitted functions do not 
have. The relative error between the the fitted functions and the results of the numerical integrations is less than 1% when the 
latter is reliable. 



In practice however the previous formulation can lead to numerical instabilities when the ratio 91/92 is large. 
Alternative forms that make use of the expansion of the logarithm when its variable is close to unity can be employed. 
They are of the form, 

0/(9i, 92) = ^oinwf^ w ( 15 (ft - 12) 4 (9i + 92) 8 (384q 4 + 2699g 2 2 g? + 384 g 4 ) B v (5, 0) 
z 1 01 841)11(7! q 2 (qi + <?2j 

-512ql& (2304 9l - UUq 2 q\ + 18103q$gf - 50908^? + 18103<$g? - 6144g|gi + 2304gf)) (63) 

and 

Pg(th,q») = 302702400^(1 +g 2 ) 4 ^ ^ ~ ^ " ^ + ^ " + 12257 ^ + U52q ^ ^ (5 ' 0) 

-512g^| (6912^ - 18432q 2 <7i + 75109g| 9l - 235924 g | 9l + 75109^ - 18432gf 9l + 6912^)) (64) 

where B v (5,0) is the incomplete Beta function of index v = j^v—vi- 
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C. Fitted forms 



As it is nearly impossible to reproduce the explicit result for the functions T and Q, we rather propose some 
analytical fits for a / and a g that are given below. 

They can be fitted with the form (where a = / or a = g), 

a a (qi, 92) = (j3 a (q 1 ,q 2 )/b a - c a (q s ) - d a (q r )) S a (q s ) (65) 

where b a is the value of j5 a for qi = q%, 

br = and b„ = , (66) 

s 429975 9 429975 v ; 

the variable q s and q r are defined by 

Qr = y/ql + q* and g s = l/"\/«r 2 + % 2 - (67) 

The functions c a (q s ), d a {q r ) and S a are fitted functions. They are given by, 

, v = 0-401228 1.46195 gs 3 

Cf[qs> 35.0987^ + 4. 13381q2 + 1 77.7967^ + 1 1 ' 

0.178036q 2 0.425873 
4.11634< ? 10 + 1 ~ 7.03063q 4 + 2.45787g 2 + 1 

0.0419954 (76.7666g 4 + 4.33008??) 0.0121732 

Sf(q s ) = ^ — + - A +0.0941118 (70) 

nqsJ 76.7666<? 4 + l 1010.66^ + 28.3738gf + 1 V ' 

0.478032^- 0.357391 

Cg ^ s ' ~ 28.4904gf + 1 ~ 42.1407^ + 1.36756^ + 1 ( > 



. . \J. -L t UUOKjq r U.^tiUUI J f69l 

f\Qr) — A 1 1CO/) 10 1 7 fl^nK!?«4 _l_ 9 AK787«2 i 1 ^ J 



. 0.348953 0.268523q 2 _^ 

9 ^ r ' ) " 2. 23826^ + l ~ 0.471174g5 + 1 { '-'' 



= 0.0351332 (9.76583^ + 3.723250 0.0580042 0.0120398 (73) 

9Vy ; 9.76583^ + 1 0.0112425g 4 + 2.06647? 2 + 1 v > 



D. The resulting propagators 



The resulting propagators are computed with the help of the Eqs. (44 45). The regular part of G 2 loop is then 
given by 



^G^ik, „)1 2 - (4vr) 2 / dft dq 2 P( qi ) P(q 2 ) «g* 



2 „2 z.2 

2 a a {q 1 /k 1 q 2 /k). 

<?1 + 92 



(74) 



It is to be noted that the integrals over qi and 52 safely converge for A-CDM type models. For a power spectrum that 
reaches a power law behavior of index n at large k one can easily check that the convergence is obtained for n < — 2 
(the condition is n < — 1 for the one-loop expression). These convergence properties will be discussed in more detail 
in the next Section. 

Finally we present the resulting shape of the propagators. At one- and two-loop level the regularized approach is 
un-ambiguous. It leads to, 



Greg— 1 — loop 
a+ 



(k,rj) = 



Regg< 1_loop 



„ + (M) exp[--k 2 a 2 d e^ 



1 



(75) 



and 



Gi e r 2 - l °°*>(k, v) = [e" + Rc <; loop (fc, V ) + ^G 2 a - lo ° P (k, r,) exp ( e 2 " 



(76) 
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FIG. 7: The resulting propagators up to two loop order in units of the linear growing mode (top) and divided by the standard 
damping factor (bottom). Calculations are made for a standard A-CDM model at z = 0. The left panels correspond to 
the density field and the right panel to the reduced velocity divergence. The (yellow) dashed lines correspond to the RPT 
expression (at one loop order) and the dot-dashed lines (green) to the 1-loop regularized expression (751, the (red) long dashed 



lines correspond to the two-loop RPT result ( 78 1 , the (blue) solid lines to the two-loop regularized expression ( 76 I and finally 



the black dotted lines to the alternative form of the two-loop regularized expression obtained from ( 43 1 



We note however that it is however possible to build a RPT like expression, e.g. to exponentiate the two-loop result 
and construct an alternative form that contains the same two-loop order terms, At 1-loop order, the RPT propagator 
reads, 



(77) 



This expression can be extended into, 



G RP T -2- 1 oop (fc; v) = eV exp f 1^2 £ 2 V + Res^-toop^ ^ 



(M) ' e^+ R ^G 2 a l l ° OP (k, V )e 



(78) 



This is a well-behaved expression since the second and third terms within the exponential is given by the second term 
of the right hand side of Eq. (58) or (59) and is negative. The resulting propagator is therefore always decreasing 
with k and with time. One can remark that if (Xf{k > 51,52) and a g (k, 91,92) vanish the expressions G^+ T (k,rj) and 
G^ T ~ 2 ~ loop (k, 77) are identical. 

These various expressions are compared on Fig. [7] It appears that in practice, for a A-CDM model the one-loop 
expressions, (77) or (75), lead to virtually indistinguishable results in the top panels. This is also to some extent the 



case for the two-loop results. Differences between the one-loop and two-loop results are however significant. We will 
in the following compare those predictions with N-body results. 

The bottom panels show more detailed comparisons of these predictions. It appears in particular that, unlike the 
RPT prediction at two-loop order, the Rog G^ + lo ° P (fc, 77) form can go negative at large enough k. This might be seen 
as a flaw (although this is of no consequences in practice because it takes place at scales that are quite damped) but 
it is corrected when one takes into account the one-loop correction to the displacement field (gray dotted lines) as it 
is the case when one uses Eq. ( 43 1 . 



15 
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FIG. 8: The regularized density propagator at 1 and 2-loop order compared to TV-body results for z = 2, 1 and 0.35. The top 
panels correspond to the propagators in units of the linear growing factor and the bottom panels to the propagators divided 
by the expected leading order damping factor. The convention for the theoretical predictions are those of Fig. [7] the results 
from TV-body simulations are shown with black points with error-bars. 



We finally give below the asymptotic forms of the multi-loop corrections in the low k limit. Formally we have 



where 



o l — loop 



G a+ (k, 77) ps e" + C i- loop fc 2 e 3 " + c 2 a - loop k 2 e 5 " 
= J dq P Q (q), 4- loop = -4tt~ J dq p o{q) , 



according to Eqs. (33) and 



2 2-loop = _ (4 ^ )2 | ^ J ^ 



More specifically, for the WMAP5 cosmological parameters we have, 



2 2 

§^ /3„(<7i,<7 2 ) Po(qi) Mto)- 



G 1+ (k,r)) w e"-" 
G 2+ (k, v ) « e*?-* 



1 



0.31 /iMpc" 1 
k 

0.17 /iMpc" 1 



^27,-2770 



o 27)-27, 



0.60 /iMpc" 1 

jfe 

0.89 /iMpc" 1 



O 4r;-4r;o 



O 4r)-4r; 



(79) 
(80) 

(81) 

(82) 
(83) 



up to two-loop order where rjo corresponds to z = 0. We note that as scales in the quasi-linear regime - say k of the 
order of O.lhMpc" 1 - the two-loop corrections induce non-negligible effects irrespectively of the shape of the damping 
tail one observes at high k. This is clearly visible on the bottom panels of Fig. [7] 



E. Comparisons with N-body results 



Finally we compare our results with N-body measurements. Here, we briefly explain the settings of the simulations 
as well as how we measured the propagators. 

We have performed iV-body simulations in periodic cubes with a side length of 2, 048/i -1 Mpc employing N — 1024 3 
particles. The initial conditions are created at z — 15 using a parallel initial condition generator developed in [241125) . 
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The code computes the displacement field based on the second-order Lagrangian perturbation theory (2LPT; [TTll26| y 
and gives the displacements as well as the initial velocities to the particles placed on a regular lattice. In doing so, we 
adopt the best-fit flat ACDM cosmology to the five-year observations of WMAP satellite [27J , and compute the linear 
matter power spectrum by CAMB [25] • We created 60 independent realizations of particle distributions, and simulate 
the gravitational evolution by a publicly available tree-PM code, Gadget2. We store the outputs at z = 3, 2, 1 and 
0.35. 

We assign the particles onto 1024 3 grid points using the Cloud-in-Cells interpolation scheme to obtain the final 
density field. This field is then Fourier transformed using the Fast Fourier transformations. The propagators are 
measured by taking the cross correlation between this field and the linear density field. Note that here we do not use 
the initial particle distribution as the linear density. Instead, we employ the seed density field which has been used 
as the source term in the 2LPT calculation. The former suffers from a slight nonlinearity, while the latter is perfectly 
linear by construction. 

We also perform simulations with different starting redshifts, box sizes, number of particles and internal parameters 
which control the force accuracy to see the convergence of the propagators. We especially pay attention to separating 
the effects of UV and IR cutoffs on the propagators, since we can cover only a finite wavenumber range in simulations 
with finite number of particles. We made sure that our box size (2, 048/i _1 Mpc) is large enough to avoid any 
systematic error larger than 1% in the propagator originating from the IR cutoff. The UV cutoff, on the other hand, 
can be a source of a small systematic decay of the propagator at small scale; this effect is almost independent on the 
output redshift, and is typically 1% at k ~ 0.3h Mpc -1 , and 2% at k ~ OAh Mpc" 1 . The starting redshift, z = 15, 
was chosen so that this systematic error is minimized given the box size and the number of particles; the force error 
(the tree force, more precisely) damps structure and leads to a stronger decay of the propagator when started at a 
higher redshift, while the 2LPT is no longer applicable at lower redshift. See Nishimichi et al. (in prep.) for more 
details. 

The propagator measured from the simulations is shown as symbols in Fig. [8j with error bars showing the 1-a 
uncertainties inferred from the variance among 60 realizations. We plot G Q +(fc) for the density field in the upper 
panels, and we also plot it after it is divided by its high-fc asymptote to make a clearer comparison. At redshifts higher 
than z — 2, the 2-loop prediction is in good agreement with the simulation data. The V-body data are slightly smaller 
than the two- loop prediction presumably due to the systematic error in the V-body simulations (see above). On the 
other hand, both one-loop and two-loop predictions can not explain the Af-body data at lower redshifts (z = 0.35 and 
and to a lesser extent z = 1). The Af-body data lie in between the two predictions. 

IV. CONVERGENCE PROPERTIES 

In this section, we are interested in the sensitivity of the results to the IR and UV modes. To this end, we introduce 
the notion of kernels that can be viewed as the response function of the propagator with respect to the initial linear 
spectrum. This notion is also engrained in the approach proposed in |29] to obtain accelerated computations of 
the non-linear power spectra. We here first properly define these quantities and explore their small and large scale 
behaviors. 



A. Monte-Carlo integrations for the three-loop contribution 

In order to gain further insights into the properties of the loop corrections, we have complemented our theoretical 
investigations with numerical evaluations of the three-loop contribution to the propagator. Unlike the two-loop 
contributions we did not try to obtain a semi-analytic expression of the results that could be used for any cosmological 
models. We rather aim here at evaluating the impact of the three-loop effects in a specific A-CDM like model. 

The results we present hereafter have been obtained from brute force calculations of G^ loop , the two-point prop- 
agator for the density starting with growing mode initial conditions. The computations focus on the calculations of 
the density late time propagator at z — 0.5 redshift. We apply the Monte-Carlo technique of quasi-random sampling 
using the library CubagS], and directly perform 2D, 5D and 8D integrations for G^ loop , G^ loop and G^ loop respec- 
tively. The transfer function of the linear power spectrum is calculated with CAMB code adopting the cosmological 
parameters determined by WMAP5 |27| . 

In Fig. [9] we separately show the one-, two-, and three-loop corrections to the two-point propagator (blue, green, 
and red, respectively). In top panel, to reduce the dynamic range, the results are divided by k 2 , i.e., G"~ loop (£;)A 2 , 
and are plotted as function of wavenumber in logarithmic scales. Note that the solid and dashed lines respectively 
indicate the positive and negative values in amplitude. The Monte Carlo results for one- and two-loop corrections 
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FIG. 9: Monte Carlo results for p (fc). The one-, two-, and three-loop corrections are respectively plotted as blue, green, 

and red lines. Top panel shows Gi+ loop (fc)/fc 2 , while the bottom panel shows the fractional difference of G"+ loop (fc) with its 
asymptotic behavior in the high k limit G"^ loop (fc) = (— k 2 <j\/2) n /n\ (indicated by black solid lines in top panel). 



accurately reproduce the analytic calculations. At low-k regime, all the plotted results show scale-independent be- 
havior, indicating that the corrections behave like G"+ loop (fc) ~ k 2 . On the other hand, at high-k regime, the three 
curves behave according to the expected forms described by Eq. ( 34 ) . They tend to approach the asymptotic form, 
G^ loop (fc) = (— k 2 (T^/2) n /nl, depicted as black solid lines. To see this clearly, the bottom panel of Fig. [9] shows 
the fractional difference between G"+ loop (/c) and its asymptotic behavior, i.e., G"+ loop (&;)/ [(— k 2 a\/2) n /n\\ — 1. At 
k > l/iMpc -1 , all the three curves consistently converge to zero. We can actually verified that the regular part of 
each contribution is finite at large k. 

In Fig. [TO] we show the relative importance of the loop corrections as a function of scale. What is plotted are the 
regular parts of the terms as defined in subsection II D It appears that at k = 0.1/i/Mpc and at z = 0.5 the one-loop 
contribution introduces a 5% correction, the two-loop contribution a further 2% contribution. At k = 0.2ft,/Mpc the 
two contributions are both at 10% level but then the mode damping is already significant. Note that as the redshift 
is increasing the importance of all these loop corrections decreases, their magnitude scaling like D 2 ^ where D + is the 
linear growth rate and p is the number of loops. 

Obviously, we hoped and expected the three-loop contribution to be smaller that the two other terms for k < l/i/Mpc 
but this is not the case. At k = 0.1/i/Mpc it gives a 10% effect! To say the least this suggests that the three-loop 
corrections, assuming they are correctly evaluated by our Monte-Carlo integrations, cannot help to match the N-body 
results at low redshift. Given the results shown on Fig. [9] we have not much reasons to doubt the validity of the 
numerical integration. To better apprehend the origin of these effects we hereafter introduce the notion of kernels. 



B. Definitions and properties of the kernel functions 

Let us define the functions K p ~ loop (k, q) as the kernels from which the contribution of the two-point propagator to 
p-loop order can be reconstructed, 

Rc <; lo ° P (fc) K^(k,q)P ( q ). (84) 
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FIG. 10: Regular parts of the density propagator Hcs Gi + lo ° P (fc) at one-, two-, and three-loop order with respectively solid, 
dashed and dotted lines. The calculations are done for z — 0.5. Note that each of this contribution scales with the redshift like 
D+(z) 2p where p is the number of loops. The light yellow regions show the parameter space where the induced corrections to 
the power spectrum are less than 1 percent. 



i^V oop (fc, q) = W ( f(q, k) + ~) (85) 



We have then for instance 

1 k 
6 q 

K^°°*(k,q) = -( 47 r)V / d ?i^2«/(f>!) P ^ ^ 

Note that the kernel functions depend themselves a priori on the initial power spectrum: K a 7 oop (k, q) is a tree order 

object, K^7 oop (k, q) a one-loop order object (and therefore a linear function of P (q)), etc. These functions give, for 
each order, the impact of a linear mode q on the amplitude of the late time mode k we are interested in. In particular 
it tells how the small-scale modes affect the large-scale modes under consideration. In the following we will focus our 
interest in understanding the high-q behavior of the kernel functions K (fc, q). 



In Fig. 11 we show the shape of the kernel functions at one, two-loop and three-loop order for k = 0.1/i/Mpc. The 



dashed line corresponds to the one-loop expression. As can be seen it is rather peaked at q rs k and we have 



Kl- lo °v(k,q)P (q) = W <Z 3 fo(?) for q « k (87) 
K 1 1 - loop (k,q)P (q) = l ^k 2 qP{q) for q » k (88) 

At two-loop order, the behaviors are qualitatively different. The function peaks rather for q = 0.5/i/Mpc, irrespectively 
of the value for k (when k < 0.5/i/Mpc). We note that 

K 2 1 ^ oop (k,q)P Q (q) ~k 2 q 2 P (q) for q » k (89) 

so that the convergence is obtained for a spectral index smaller than —2. This corresponds to the result mentioned 
in the beginning of subsection |III D| These trends are amplified for the three-loop results shown with a dot-dashed 
line for which an even lower power law index is required for convergence. In general the convergence properties of the 
multi-loop kernel are determined by the properties of the functions F n {(\i) and G n (qi) and how they behave when 
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FIG. 11: The shape of the kernel functions P (q)K 1 - loop (k, q) (blue solid line), P (q)K 2 - loop (k, q) (green dashed line) for 
k = 0.1/i/Mpc and P {q)K 3 ~ loop (k, q) (red dotted line) as a function of q for z = 0.5. 



one or their argument is, in norm, much larger than the sum of the wave modes. As mentioned in [30] it is to be 
noted that the Galilean invariance of the motion equation implies that, 

F n (qi,...,q„) K^— when q, > |$^<fe|, (90) 

Ql j 

whenever one of the is much larger than the sum. This can be seen at an elementary level on the properties of the 
vertex function a(ki,k2) and /3(ki,k2): they both vanish when the sum of the argument goes to 0. The property 
(901 has direct consequences on the properties of the loop corrections. As a result, the loop correction takes indeed 
the form, 

GP- lo °P(£:, rf) = (4^ j qfdq, . . . q 2 p dq p P (<Zi) • • ■ Po(q P ) -F p - loop (fc, Qu---, %) (91) 

with 

J-P- lo °P(fc, qi , . . . , q p ) = ^^y' J d 2 <Ji • . . d 2 g p F 2p+ i(k, qi) - qi , . . . , q p , -q p ) (92) 

where the integration is here restricted to the wave mode angles. The property ( 90 1 then implies that 

k 2 

T p - loop (k, q u . . . , q p ) ~ -g when q t > k. (93) 

A simple and effective way to capture this property is to write F(k, q±, . . . , q p ) following the functional form, 

fc 2 

F{k, qi, . . . , q p ) = =; — ^a p (k/qi, k/q p ) (94) 

generalizing Eqs. 1 58|59 1 and where a(k/qi, . . . ,k/q p ) is always finite. This property determines the converging 
properties of the p— loop correction to the propagator. Clearly the converging properties of the loop corrections 
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deteriorate. It can easily be shown that formally, for a power law spectrum of index n s , the convergence is ensured 
only when, 

2 

n s <-3+-. (95) 
P 

It clearly appears that the convergence domains shrink with the order. It leads to a greater sensitivity to the large 
wave-modes as one increases the loop order. What we stress here is that the kernel functions can explicitly be 
computed and thus one can exactly quantify the impact of large wave modes to each p-loop correction terms. 
Note that it is possible to define in general the kernel of contributing terms as 

= / J Kt b (k,q)P (q) (96) 

where (J denotes a peculiar term. One can then compare the kernels of the different terms that contribute to the 
non-linear power spectrum. It is worth mentioning that the kernels for the two-point propagators are the broadest. 



Indeed any other kernel function in the large-g limit is going to take advantage of the asymptotic form ( 90 ) at least 
two times instead of one so that 

Kl + (k,q)P Q (q)~^ q 3l P (q) for q » k (97) 



where I is the number of loops appearing in the diagram jj. This is to be compared with Eqs. (88 89). 

Clearly the loop corrections to the propagators are then the contributions built from perturbation theory calculations 
that are most sensitive to the small scale modes. At low redshift our findings imply that a regularization procedure in 
the UV domain should be used in order to make the contributions of these terms realistic. Unless partial resummation 
of higher order terms could provide us with a self-consistent regularization scheme, it should be derived from external 
pieces of information. This regularization could be associated for instance from the development multi-flow regimes. 
For instance it was advocated in |16j that the impact of the small scale physics could be captured with an effective 
theory approach where effective anisotropic pressure terms are introduced in the fluid. Our results however suggest 
that such effects are subdominant at the one-loop order - the contributions of such terms need not be recalculated 
because they are naturally protected enough from UV physics - but not at two-loop order for z < 0.5. It also suggests 
that when z gets larger the importance of such corrections should vanish away even at 2 loop order at the nonlinear 
regime starts at larger wave-modes. 



V. CONCLUSIONS 



Two-point propagators are key ingredients for the computation of power spectra in a cosmological context. These 
are also quantities for which it is easier to make analytical predictions with a Perturbation Theory approach in parts 
because they are expected to be damped in the large-wave modes limit in parts because the number of terms in the 
loop corrections are more easily computable. In this paper we have presented the most advanced results describing 
the expected behavior of the two-point propagators from a Perturbation Theory point of view. In particular we have 
exploited and extented the forms presented in [T] to propose theoretical expressions for the two-point propagators 
that incorporate up to two-loop results. The examination of these forms led to the concept of regular parts for loop 
contributions that are found to be finite even when the r.m.s. of the displacement field, <7d, goes to infinity: regular 
parts are found to be the parts that contribute to the expression of the propagators once the dependence with a& has 
been explicitly taken into account. They also can be used as a consistency test for Monte-Carlo numerical calculations 
of the propagators as they predict the leading k behavior of the contributing terms. 

The main part of this paper is the explicit calculation of the two-loop contributions to the two-point propagator. 
We were not able to produce a simple analytical form but we rather present, in section |TiT| simple fitting functions 
that capture all the expected properties of these contributions, in particular their asymptotic properties. These forms 
can be used for any cosmological models. Note that we propose different variants for the two-loop expression of the 
two-point propagator but they differ only mildly. There results are presented in Fig. [7] This is the main result of 
the paper. It has been compared with iV-body results. They show that two-loop corrections help significantly in 
predicting the shape of the two-point density propagators at z > 0.5. 

For lower redshifts the two-loop terms appear to be too large when compared to iV-body codes. We found that a 
possible reason for such a discrepancy is due to the increasing sensitivity of the loop corrections to the small-scale 
modes. A significant fraction of the contributions indeed originate from modes that are in the fully nonlinear regime. 
At those scales the whole perturbation theory scheme breaks down as the fluid has not only evolved into the nonlinear 
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regime but also contains multi flow regions with order unity anisotropic pressure. The latter is not captured by the 
motion equations we use and should induce extra contributions as it has been put forward in recent papers, p~5|. I16j. 
We do not think however that it has significant impact on the one-loop results. Those effects are rather expected to 
be significant at two-loop order and for z < 0.5. This sensitivity to the small scale physics is shown to be stronger for 
the three-loop contribution that we computed with the help of Monte-Carlo integrations. 

About the converging properties of the Perturbation Theory series we note that loop corrections, at any order, 
induce large-scale k 2 corrections to the propagators and consequently to the power spectra. They are explicitly given 
in Eqs. ( 79|83 ) up to two-loop order with coefficients that depends on the overall shape on the linear power spectrum. 
Again it stresses the relative importance of the sensitivity of the standard Perturbation theory calculation to the UV 
domain. We also note that this sensitivity is controlled by the Galilean invariance which leads to the general property 
given in Eq. (90). This form appears to be a key ingredient for the converging properties of the Perturbation Theory 
series. It indeed sets the converging properties of both the diagrams that correspond to propagator loop corrections, 
as given by Eqs. (88 89), and the other ones, as depicted in Eq. (97). Clearly if this property is dropped, the 



sensitivity of loop corrections to the UV domain is more severe. We note that if in standard gravity cosmologies, 
the relation (90) is satisfied, this is not necessarily so in cosmological models incorporating modified gravity effects 
or in general a clustering quintessence component. This is the case for instance in models explored in 31-33 that 
incorporate extra dynamical degrees of freedom |41j . 
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